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We address the question of how a quantum computer can be used to simulate experiments on 
quantum systems in thermal equilibrium. We present two approaches for the preparation of the 
equilibrium state on a quantum computer. For both approaches, we show that the output state of 
^\ . the algorithm, after long enough time, is the desired equilibrium. We present a numerical analysis 

of one of these approaches for small systems. We show how equilibrium (time)-correlation functions 
On : can be efficiently estimated on a quantum computer, given a preparation of the equilibrium state. 

The quantum algorithms that we present are hard to simulate on a classical computer. This indicates 
that they could provide an exponential speedup over what can be achieved with a classical device. 
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Cn ' I- THE LIMITS OF CLASSICAL COMPUTATION 

> 

The power of quantum computers has been demonstrated in several algorithms, of which the most striking have 
been Shor's factoring algorithm and Grover's search algorithm |||. From the very start however, the quantum 
computer has also held the promise of being a simulator of physical systems. This is the content of the physical version 
of the Church- Turing principle proposed by Deutsch Q . Thus we might expect that the universal quantum computer 
can be used to simulate any experiment that we could do on a real physical system. If such a simulation can be done 
efficiently (that is, without exponential slowdown), it is clear that this could be one of the major applications of a 
quantum computer. This promise seems to have been only partly fulfilled until now; it has been shown by several 
researchers |||| that a simulation of the unitary time evolution of a physical system that possesses some degree of 
! 1 * . locality (which realistic physical systems do) can be accomplished efficiently on a quantum computer. However, many 
quantities of interest that are determined by experiment, or by the use of classical simulation techniques, relate to 
open quantum systems, in particular to systems in thermal equilibrium. The thermal equilibrium (Gibbs) state (in 
the canonical ensemble) of a Hamiltonian H is given by 



N 

W = ^i— |m)H, (1.1) 
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where \m) (E m ) are the eigenvectors (eigenvalues) of H. Z is the partition function 



N 

Z (1.2) 

and (3 = ™ where k is Boltzmann's constant and T the temperature. The physical systems that concern us in this 
paper will have a finite dimensional Hilbcrt space Ti that can be decomposed as 

H = Hi®H 2 ® ...®H n , (1.3) 

where each Tii represents a small, constant Hilbert space, typically associated with some (generalized) spin or other 
local degree of freedom. The Hamiltonian couples these local Hilbert spaces, for example in correspondence with a 
rf-dimensional spatial lattice, so that there is only coupling between adjacent "spins" on this lattice. The quantities 
of interest, computed in experiment or in a classical computation, are of the form 

Tr Oi(ti)0 2 (t2)0 3 (i3) ■ ■ ■ O k (t k ) Pp , (1.4) 

where 0j(ij) are (possibly time-dependent) observables. Both for classical systems as well as for quantu m sy stems, 
computational Monte Carlo methods have been developed to estimate correlation functions as in Eq. (1.4) 0-^). 
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The quantum Monte Carlo method for systems at finite temperature relies on a transformation introduced by Suzuki 
ijiof that maps an initial quantum system on a cZ-dimensional lattice onto a (d + l)-dimensional classsical system. 
This conversion then ma kes it possible to use classical computational sampling techniques to estimate correlation 
functions as in Eq. (|l.4|) . There seem to be (at least) two situations when this approach runs into trouble and no 
good computational alternatives are available 01: (1) the correlation functions depend explicitly on time t, and (2) 
the quantum system is of a fcrmionic nature. We will give a short explanation of why these problems arc encountered. 

The transformation from a classical to a quantum system is based on the generalized Trotter formula. Let H = 
53j=i Hi where each Hi is a Hamiltonian on a small constant Hilbert space. The Trotter formula reads 



e 



aH = lim ( e °Hi/n e aH 2 /n e *H k /n\ ^ 



The partition function Eq. ( |l.2| ) (and similarly correlation functions as in Eq. ( |l.4| )) can be rewritten, using the 
Trotter formula and the identity J2 m l a i")) ( a i"? I = 1j where the pair of indices (i, j) labels a choice of basis, as 

Z = Tre-P B = ]T P{ai:j} , (1.6) 
where p^ a . .\ is a distribution over the values of the collection of variables {flij} and j indexes the repetitions of the 



factors of Eq. (1.5) from 1 to n. If the distribution is nonnegative, we can write p{ a . .y — e He ff^ ai ^ where -£f e // is 



now a classical Hamiltonian given by 

n k 



H eff({ai,j}) = n 1 ™ o XlX)^ i ( ai J' ai + 1 'j)- ( L7 ) 

with afc+ij = aij+i, ak+i, n = ai,i and 

Hi(a,b) = ]og((a\exp(-PHi/n)\b)). (1.8) 

The distribution py a . y will only be nonnegative when the matrix elements (a\ exp(—(3Hi/n)\b) are positive. Thus 
it is important to choose the right sets of basis states |a™) to make the conversion to a classical sampling problem 
with a positive distribution. There are fermionic systems such as certain Hubbard models |9| in which it does not 
seem to be possible to choose such a good basis. For these systems it has turned out to be very hard to get good 
estimates of correlation functions by using classical Monte Carlo techniques. This problem is usually referred to as 
the "sign" problem. 

When we are to compute time-dependent quantities, for example the fun ction f(it) = Tr e tHt Oie~ lHt O2P0, we 



need to use an imaginary time r = it to perform the conversion of Eq. (L5) to a classical system (we expand e lHt 
with the Trotter formula). From the classical Monte Carlo sampling of the function /(r) for real r, we estimate f(r) 
and then we could in principle analytically continue this function. However, we only have a finite number of samples 
of the function and each sample point has some inaccuracy. The errors that are introduced in estimating the Fourier 
components f(to) from this data give rise to large fluctuations when we reconstruct f(it) with the Laplace transform 

f(it)= f°° dwe-^/H, (1.9) 

J — oo 

resulting in a bad approximation for the time correlation function f(it). 



The relevance of estimating a simple time correlation function (an example of Eq.(L4)) such as 

Tr [A(t), B(t')] Pp = ([A(t),B(t')]) s , (1.10) 
where A and B are some Hermitian Heisenberg operators of the system, cannot be overestimated. Let us recall the 



many contexts in which Eq.(l.lO) is used in describing experimental properties of many-particle quantum systems 

When A = B = u, where u is the displacement field of a crystal, fll.lOD describes the phonon dynamics of solids as 
probed by inelastic neutron scattering. When A and B are the number-density operator, the dielectric susceptibility 
is represented; this correlation function describes a variety of other experiments, including x-ray photoemission and 
the so-called x-ray edge singularity. When we study the current-current response function, we obtain the electrical 
conductivity as described by the Kubo formula. (The density-density and current-current response functions are 
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intimately related via the continuity equation.) Spin-dependent quantities are also of interest: with the spin-spin 
correlation function, information is obtained about the magnetic susceptibility, and thus the magnon dynamics of 
ferromagnets and antiferromagnets, the Kondo effect, and the magnetic-dipole channel in neutron scattering. And 
finally, if A and B involve anomalous pair amplitudes which involve Fcrmion operators like aj.(fc)of(— k), the presence 
and dynamics of a superconducting phase can be probed. 

In short, the dynamic pair correlation functions provide a window on many of the interesting quantities in exper- 
imental physics, and it would be highly desirable to have a method of obtaining estimates for these quantities by 
simulation on a quantum computer. We will present some methods below for doing this. 

In this paper we develop an approach to tackle these problems on a quantum computer. We break the problem 
into two parts: First, we present an approach to prepare our quantum computer in the equilibrium state pp of a 
given Hamiltonian (sections || and III). We will give two alternative routes to prepare an equilibrium state. Next 
we describe a procedure for efficiently estimating quantities as in Eq. ( jO| ) given that the equilibrium state has been 
prepared (section IV). We will not attempt to prove that our algorithms run in polynomial time even for a certain class 
of quantum systems H and/or for certain ranges of (3. The equilibration problem, in its full generality, is expected 
to be a hard problem. Even classically there is a large class of systems that exhibit a feature called frustration, for 



which calculating the partition function Z as in Eq. (1.5) is a P'-complete problem [Q. Also, for these systems, 
deciding whether the energy of the ground state is lower than some constant K is an NP-complete problem JjjJ . The 
quantum problem has an added difficulty: We cannot assume that we know the eigenvectors (and eigenvalues) of the 
Hamiltonian of the system that we would like to equilibrate. There has been no demonstration yet that a quantum 
computer can exponentially outperform a classical computer in estimating the partition function for certain classical 
systems, which would enable us to sample efficiently from the classical Gibbs distribution p^). 

The quantum algorithms that we present are hard to simulate on a classical computer. In both of our equilibration 
algorithms we use the fact that one can implement the unitary time evolution of a local Hamiltonian on n qubits in 
a polynomial number of steps in n on a quantum computer fa] . A direct simulation of this procedure on a classical 



computer would cost exponential (in n) space and time and is therefore unrealistic. As we will show in section |IV| , 
given a preparation of an equilibrium state, there exists an efficient procedure on a quantum computer to calculate 
(time-dependent) correlation functions. As we discussed above, there is no general efficient classical algorithm with 
which one can estimate time-dependent correlation functions. Our quantum algorithm provides such an algorithm 
for a quantum computer. Lloyd and Abrams ]l6| ] have shown that the unitary simulation of a fermionic system such 
as the Hubbard model, either in first or second quantization, can be performed efficiently on a quantum computer. 
The quantum algorithms that we will present will use this unitary evolution as a building block. Therefore these 
algorithms can be used to compute correlation functions for the Hubbard model on a quantum computer. This is a 
task for which we do not have a good classical algorithm, due to the "sign" problem, as we pointed out above. 

We focus our efforts on quantum equilibration algorithms for Hamiltonians of which the eigenvalues and eigenvectors 
are not known beforehand. These are the Hamiltonians of, for example, Heisenberg models (in more than two 
dimensions), Hubbard models, t-J models, XY models, or many-electron Hamiltonians in quantum chemistry. On the 
other hand, knowing the eigenvectors and eigenvalues of a Hamiltonian, such as in the Ising model, is no guarantee 
that there exists an efficient (polynomial time) classical algorithm that produces the equilibrium distribution. The 
situation is simil ar fo r quantum algorithms; we do not know in what cases the equilibration algorithms presented 
in section |n] and III give rise to a polynomial time algorithm (see also [ jl5[ for quantum algorithms for Ising-type 
models). 

The process of equilibration is also essential in the actual realization of a quantum computer. One of the assumptions 
underlying the construction of a quantum computer |jl7f is the ability to put a physical system initially into a known 
state (or a thermal equilibrium state in the NMR quantum computer the computational 1 00 . . . 0)(00 . . . 0| state. 

The way this is done in an experimental setup is to let this state be the ground state of a natural Hamiltonian and 
subsequently to cool to low temperature such that the probability of being in this ground state is some constant. 
This natural Hamiltonian must be sufficiently simple for this equilibration to be achievable efficiently and also be 
sufficiently weak or tunable not to disturb the computation later on. 

II. EQUILIBRATION I 
A. Introduction 

The canonical ensemble is the ensemble of states {pi, or a density matrix p = X^PilV'iKV'il! such that p has 

a given energy-expectation value 

Tr Hp = {E). (2.1) 
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The equilibrium state in this ensemble (Eq.( |l.l| )) can be obtained by maximizing the von Neumann entropy of p 
under this energy constraint. Another way in which the canonical ensemble is defined is by considering the possible 
states of a system that is in contact with an infinite heat bath at a certain temperature T. The total energy of system 
and bath is constant, but bath and system exchange energy, so that the system equilibrates. This directly suggests 
that the way to prepare the equilibrium state on a quantum computer is to mimic this process. In considering the 
computational complexity of such a procedure, we will have to include the space and time cost of the bath, which may 
be large. Also, the intuitive picture of equilibration between a weakly coupled large bath and system does not tell 
us anything about the rate at which this equilibration occurs. Furthermore, the equilibration process assumes a bath 
that is already in its equilibrium state. Can we make the bath simple enough that this bath state can be prepared 
efficiently? In this section we study this process of equilibration. We present an algorithm and we derive expressions 
that completely characterize the equilibration process in an idealized case: the coupling between the bath and the 
system is very small, the bath is very large, and the time of interaction is large. We then proceed by a numerical study 
of the algorithm in realistic cases where the bath is of finite dimension, the strength of the interaction is non-zero, 
and the interaction time is limited. 



B. The algorithm 

Definition 1 Equilibration algorithm I. 
Input-parameters: 

-H s , the Hamiltonian of a N = 2 n -dimensional quantum system. 
-(3 , the inverse temperature. 

-Hb, the Hamiltonian of a K = 2 k - dimensional "bath" quantum system. 

-\H s b, where H s b is the N K -dimensional "bath-system" interaction Hamiltonian and A is the parameter that measures 

the strength of the interaction between bath and system. 

-t, the interaction time between bath and system. 

-r, the number of times the bath is refreshed in the algorithm. 

Define the total Hamiltonian of system and bath as 

H = H s ®l K + l N ®H h + \H sh , (2.2) 

and the trace-preserving completely positive (TCP ) map S\,t os 

SxAp) = Tr b e iHt p®p b ,(3e- im . (2.3) 

1. Prepare system. We prepare then qubits in the computational state: 1 000 . . . 00) (000 . . . 00 1 . 

2. Prepare bath. We prepare the k qubits of the bath in their equilibrium state pb y p of Hi,. 

3. Evolve system and bath for time t and discard bath, that is, perform the superoperator S\,t of Eq. /l2.t\) . 
4- Repeat steps 2 and 3 r times such that 

|| S^flOOO . . . 00)(000 . . . 00|) - Sj^flOOO . . . 00)(000 . . . 00|) \\ tr < e, (2.4) 

for all r > r and e is some accuracy. See Appendix p| for the definition of \\ . \\ tr . 



We put several constraints on H Sl Hb, and H s b- We will use local Hilbert spaces as in Eq. (1.3) of dimension 2 
(qubits). H s must be a "local" Hamiltonian. We define a c-local Hamiltonian on n qubits as one that can be expressed 
as 

poly(n) 

H s = 1 N/c®h l , (2.5) 

i=l 

where each hi operates on a tensor product of several small qubit Hilbert spaces, whose total dimension is c. We 
will also assume that the eigenvalues of H s are all distinct; the spectrum is non-degenerate. This will simplify the 
upcoming analysis. In order to treat Hamiltonians with degenerate spectra a change in the perturbation theory of 



Section [ID will have to be made. We expect however that with that change the main result of Section HE, namely 



succesful equilibration in the idealized case, will still hold. H s b has the linear coupling form 
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S®B, 



(2.6) 



where both S € B(TC S ) and B e B(H.b) are local Hamiltonians. Hi, is the Hamiltonian of a system of non-interacting 
qubits, i.e., it is a sum of single-qubit Hamiltonians: 



if/2 



(2.7) 



The bath's equilibrium state factorizes into a tensor product of qubit equilibrium states associated with each hi. 



(2.8) 



This enables us to prepare the bath (step 2) efficiently. Appendix |B] shows that it will cost 2k elementary qubit 
operations to perform step 2. The locality of H s , Hb, and H s b is required in order to be able to simulate the unitary 
time evolution e tHt in time proportional to t 2 /S where S is the accuracy with which gates are implemented |5 19 1. 
We also choose 



(B) b = Tr Bpb.p = 0. 
To understand the effect of a non-zero (B)b we rewrite H as 

H = (H s + X(B)bS) ®l K + lN®Hb + \S® B', 



(2.9) 



(2.10) 



where (B')b = 0. Thus choosing a non-zero (B)b effectively corresponds to a change in the Hamiltonian of the system. 
We now discuss the last step of the algorithm, step 4. When the superoperator Sx,t has the equilibrium state p s ^ as 
its unique fixed point, then Eq. (2.4) for all r > tq implies that that 



ST t (|000 ... 00) (000 . . .00|) - p, tP \\ tr < e. 



(2.11) 



for all r > ro, that is, the equilibration process leads to succesful convergence to the equilibrium state. There does 
however not exist a straightforward implementation of step 4. The first problem is that we would have to check 
the closeness of the rth and the (r + l)th iteration of S\j for all r > tq. In practice this has to be replaced with 
choosing a finite set of iterations r for which the invariance of S r (|00 . . . 0) (00 . . . 0|) is tested. This problem is also 
encountered in classical Monte Carlo simulations. The second problem, which is a purely quantum phenomenon, is 
that by measuring p r = t (p) we might disturb p r . Thus to compare p r with p r +i we would have to run S again 
for r + 1 times. To assemble some statistics on the difference between p r and p r +i we have to run r iterations of S 
several times. These considerations about the verification of the convergence of the equilibration process are of course 
not special to the use of a quantum computer; they are the same as in the equilibration of a quantum physical system 
in an experimental setup. Furthermore, it would be an impractical task to try to measure all the matrix elements of 
p r ] p r contains an exponential amount of data of which we can extract only a polynomial amount by measurement in 
polynomial time. The best way to proceed is the same as what one does in classical Monte Carlo simulations ||. If 
the goal of the computation is to estimate Tr Op s ^ then one assembles the datapoints 



O r = Tr Op r 



(2.12) 



until \O r — O r +i\ < e for a sufficiently large set of iterations r > tq. The same procedure can be carried out when 
the goal of the equilibration is to compute a time-dependent correlation func tion such as Eq. (1.4). 

In the remainder of this section we will analyse this algorithm. In section [I C we give some general properties of 
TCP maps. In section II D we disc uss th e non-hermitian perturbation theory that will be the basis of the analysis of 
S\ y t in the idealized case. In section HE we derive explicit expressions for the idealized case. The idealized case is the 
case obtained by taking the li mit s A — > 0, k — > oo and t — > oo. We develop a perturbation theory on the basis of the 
assumption that S\j of Eq. (2.3) is diagonalizable. Then we can sh ow t hat i n th is idealized case the process has a 
unique fixed point which is the equilibrium state. Finally, in sections [I G and II H we present results from numerical 



simulations in realistic cases. The following questions will be adressed: 

1. How does k, the number of bath qubits depend on n, the number of system qubits? Are they polynomially related? 

2. What is the influence of different choices for H b , S and B (Eqs.(2.6),(2.7))? 

3. How do the r, A, and t required for successful equilibration depend on n generically? 

The dynamics of open quantum systems, like the system in our algorithm that interacts with a bath, are most often 
studied with the use of a generalized m aste r equation [ p0[ . The exact master equation in integral form describes the 
time evolution of p(t) = S\ :t (p) of Eq. (2.3): 
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p(t) = e~ iCst p{0) - A 2 / df f dt"e~ lC ° {t ~ t ' ) M{t' 1 t")p{t"). (2.13) 
Jo Jo 

where C, the Liouvillian, is defined as 

C(p) = [H,p}. (2.14) 

so that C s (p) = [H a ,p] etc. The operator Ai(t',t") is the "memory kernel", 

M(t',t") = T^C^c-^-^^'-^C^pk. (2.15) 

The form in which the master equation is most often used, however, is one in which two simplifying approximations 
are made: (1) the Born approximation. This relates to the weakness of the interaction parameter A. (2) the Markov 
approximation. The process described by S\j is Markovian if we can write 

S x .t{SxAp))=Sx.t+s(p). (2.16) 

This typically occurs when the rate at which the effect of the system on the bath is erased in the bath (in the sense of 
being spread throughout the bath) is much faster than the rate at which the system evolves; the system sees a "fresh" 
bath every time. In our algorithm this loss of correlations in the bath is enforced when after some time t the bath is 
replaced by a new bath (step 4). We would not be able to truly equilibrate a finite system with a finite-dimensional 
bath if we would not keep refreshing it. Since the global evolution of bath and system is unitary, eventually we will 
get back to the initial unentangled state and, after tracing over the bath, to the initial state of the system (a so-called 
Poincare recurrence). Whether Markovian dynamics is justified will depend on the size of the bath, the strength of the 
interaction and the length of the interaction time. There are ways to make a simple but naive Markov approximation 



in Eq. (2.13) that lead to a master equation that fails to describe TCP dynamics [|2lJ|22j. The form of the master 
equation that does incorporate both the approximations and yields a physical completely positive map is the master 
equation in Lindblad form p3| : 

^ = -i[H s ,p(t)]+Lp(t) (2.17) 
where L [|||22 can be expressed with a basis of operators Fi as 

N 2 -l 

L P (t) = - aki([F k p(t),F?} + [F k ,p(t)F?]), (2.18) 

fc,2 = l 

where <z&; is a positive semi-definite matrix. In a Lindblad equation describing the equilibration process, we expect 
L to depend on the system Hamiltonian H s . The equilibrium state p s ^ - if the algorithm is successful- should be a 
stationary state of the process, which implies that [H s ,p s jj] = and 

L Ps ,p = 0. (2.19) 

Davies [pi|-[26f has demonstrated that a process described by Sxj where the bath is an infinite-dimensional quantum 
system (for example a quantum field) does equilibrate any quantum system in the limit where A — *■ 0, t — > oo, but X 2 t 
stays constant. By carefully taking a Born and Markov approximation, he derives a Lindblad equation of the form 



such that Eq. (2.19) is obeyed. We will perform a similar analysis here. The main point of difference is that we use 
a perturbative analysis of the dynamics which is only valid for small A 2 t, but coincides in this regime with Davies' 
result. We furthermore obtain more explicit expressions for the dynamics in this limit. 

One can write the most general form o f an L that obeys a quantum detailed balance p?j condition, a stronger 



requirement that the stationarity of Eq. ( 2.19 ). Now, one might ask the following question: Could we implement 
this corresponding superoperator directly, without the use of a weakly coupled large bath, so as to save us time and 
space? We believe the answer is no, as L will depend on the eigenvectors and eigenvalues of H s , which we do not 
know beforehand. 



C. Some useful properties of TCP maps 



In this section, we study some essential properties of the superoperator defined as in Eq. (2.3). This superop- 
erator is a TCP map 



G 



Sx.t : B(H N ) B(H N ), (2.20) 

where B is the algebra of bounded operators on the Hilbert space Hn- The set TCP[iV, N] is the set of TCP maps 
S: B(H N ) -> B(H N ). 

The elements of B(Hn) can be represented as N x N matrices. An alternative and convenient way to represent 
B(Hn) is as a iV 2 -dimensional complex vector space C N 

I : X & B(H N ) ^ (xh £ C N * . (2.21) 

This representation leads to a matrix representation of a TCP map S on C n2 . Let Ai be the operation elements of 
S, i.e. 

S( X )=Yl A iXAl Yl A Ui = lN. (2.22) 

i i 

Then 

x'mn = (<S(X))mn =^^(A)mk{x)kl(A\)ln = ^S mnM {x)kU (2.23) 
i k,l k,l 

with 

Smn.ki = ^2(A)mk(A^)i n . (2.24) 

i 

One can then study the eigenvectors and eigenvalues of the matrix representation of a TCP map. First, we will 
give three useful properties of TCP maps that follow directly from their definition: 

Property 1 Let B pos G B be the set of positive semi- definite matrices. Let S G TCP[iV, N]. Then 

p G Bpos =► S{p) G B pos , (2.25) 

as S is (completely) positive. Let x be an eigenvector of S with eigenvalue //, S(x) = PX- We have 

Trx^0=^ n = l, (2.26) 



as S is trace-preserving. Let Ai be the operation elements in the decomposition of S as in Eq. (2.2% ). If x * s an 



eigenvector of S with eigenvalue p, then x^ is also an eigenvector of S with eigenvalue /i*. This follows from 

(<S(x))t=£(A iX 4)t=«S( x t). (2.27) 



Let B pos i be the set of positive semi-definite matrices that have trace 1, i.e. the density matrices. Thus Property 
|l| implies that if a density matrix p is an eigenvector of the superoperator, it must have eigenvalue 1, that is, it is a 
fixed point of the map. On the basis of the TCP property of a map S, we can also show the following 

Proposition 1 Let S € TCP[iV, N}. All eigenvalues p, of S have \p\ < 1. 

Proof (by contradiction): Assume x is an eigenvector of S with eigenvalue > 1. Note that Property |l| implies 
that x h as Tr x — 0- If X is hermitian, p will be real. As x is traceless, it must have at least one negative eigenvalue. 
One can always find a density matrix p and a small enough e such that p' = p + ex is still a density matrix. Let S 
operate r times on this density matrix. For large enough r the result S r (p + ex) = S r (p) + e l J - r X wml n0 longer be 
a positive semi-definite matrix: take the eigenvector \tp) of x corresponding to the lowest (negative) eigenvalue A m i n . 
Then 

(i>\S r (p)\i>) + ep r (^\ X m < 1 + ep r X min , (2.28) 

will become negative for large enough r. But Property |l| implies that S r {p') is a density matrix, thus \p\ cannot be 
larger than 1. When x is non-hermitian, we reason similarly. One can find a density matrix p and a small enough e 
such that p' = p + e(x + X*) is a density matrix. Let S{x) = PX = ImI^X- L et A m i n , r be the smallest (and negative) 
eigenvalue of the traceless hermitian matrix e l ^ r x + e~ l ^ r x^ ■ Then 

MSVM = (MS r (pM + eWme^X + e-^M < 1 + e | ^ T A min , r , (2.29) 

will become negative for some large r (A m i n . r is a quasi-periodic function of r so it cannot be small for all large r). □ 
Another property about the existence of fixed points can be derived: 
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Proposition 2 Let S G TCP[AT, N]. S has a fixed point (which is a density matrix). 

Proof. The set of density matrices -B pos ,i G B(Hn) is convex and compact. S is a linear continuous map and 
S(p G -B P os,i) G B P os,i- Then the Markov-Kakutani Theorem V.10.6 of |22| applies. □ 

The existence of a fixed point does not by itself guarantee that the process described by S is "relaxing" , that is 
linv^oo S r (p) — po for all p where po is the fixed point. The existence of such a limit depends on whether the fixed 
point is unique. This following Proposition proves that when there is unique fixed point, relaxation will occur and 
the relaxation rate is determined by the second largest eigenvalue of S |29[ : 

Proposition 3 Let po G B pos ^(H.n) be the unique fixed point of a TCP map S. Let \k\ = max m | AIm ^ 1 \p m \, the 
absolute value of the second largest eigenvalue ofS. Then for all density matrices p we have 

II S r (p) - po ||tr< CWpoly(r)|K| r . (2.30) 

where C'n is a constant depending on the dimension N of the system and poly(r) denotes some polynomial in r. Thus 
for all density matrices p 

lim || S r (p)- Po \\tr=0. (2.31) 



Proof. Let pi be the eigenvalues of S. Let s be the number of distinct eigenvalues. We can bring any matrix S into 
Jordan form J by a similarity transformation M pl|: 



S = MJM~ 



where 



J 



i=l 



QnPi + Ni). 



(2.32) 



(2.33) 



Pi are orthogonal projectors and Ni is a matrix of Is above the diagonal in the ith block or Ni is the matrix. When 
the eigenvalue pi is nondegenerate Ni is the matrix. We therefore have NiNj — for i ^ j an d PiNj = for i ^ j. 
Call the unique largest eigenvalue po = 1 and the corresponding projection Pq. As in Eq. (2.32) one can write 



S r = MJ r M~ 



(2.34) 



where J r equals 



(2.35) 



where N- is a nilpotcnt matrix in the ith block whose matrix elements are all smaller than or equal to rp\ . Note that 
A^o is not present as p is unique. Let 5° be MPqM^ 1 or S°(p) = p . We use || A||t r < V~N \\ A\\ 2 . Note that || A\\ 2 
refers to the Euclidean norm of A represented as a vector. This follows from (5^i=i l^il) 2 — \ x i\ 2 ^ or complex 

numbers Xi. We have first of all 



\\S r (p)- Po \W<VN\\ (S r -S°)(p) || 



(2.36) 



This expression can be bounded with the use of the similarity transformation M to 

II S r (p) - po \\ tr < VN \\\M(S r - S^M- 1 |||a < C hN \\\J r - P ||| 2 (2.37) 

whereJJL |||2 is defined in Appendix |a| and we use || p 1 1 2 = Tr p 2 < 1 for density matrices. Using the expression for J r , 
Eq. (2.35), we can also bound 



\J r -Po His < poly(r)C 2iJV |K| r 



(2.38) 



Combining Eq. (fTf?]) and Eq. gives us the desired result Eq. ( |2 . 3C|) . E q. ( |2.3l| ) then follows as \n\ < 1 by 

Proposition |l|. If S is diagonalizabl e, th e nilpotents Ni in expression Eq. (2.35) are not present. By going through 
the same steps, a bound as in Eq. (2.30) can be derived without the factor poly(r). □ 
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We refer the reader to p2] for discussions and references concerning the existence of a u nique fixed point and other 
properties of relaxation for a process that is descr ibed by a Lindblad equation, Eq. ( 2.17 ). 

The bound on the rate of convergence of Eq. (2.3C) is far from optimal for small r as we know that for any two 
density matrices p\ and pi, \\ pi — P2 ||tr< 2. However, it is not so bad as to invalidate the main conclusion that 
one would like to draw from it. If S is diagonalizable and |k| = 1 — a/n c then a (polynomial) number of iterations 
r = — (In 1/e + In Cat), for large n, results in 

\\S r (p)- Po \\ tr <e, (2.39) 

where we used lim m ^ 00 (l — x/m) m = e~ x . If S is not diagonalizable the convergence is possibly slowed by the factor 
poly(r), but there still will be a polynomial relation between |k| and r. 

Finally we give a result which relates members of TCP[iV, N] to the stochastic matrices. A real matrix M is 
stochastic when the entries of its columns add up to 1, i.e. J^- Mij = 1. 

Proposition 4 Let S £ TCP[iV, N]. 5 mmjlm G R, and, Vn, ^2 m S mm . nn — 1; that is, the elements S mm . nn form an 
N x N stochastic matrix in the indices m and n. Also, VVt, k, n 7^ k, *^2 m S mmyn k = 0. 



Proof: S mm ,nn <= P follows directly from Eq. ( |2.24| ). For the rest, we impose the unit trace condition on Eq. ( 2.23 ): 

1 = S mm ,kiPki- (2.40) 

This must be true for all density matrices represented by p. Taking pki = <5fc,z<5fc,fc gives the desired result 

1 = ^ Smm,k k - (2-41) 
m 

We now separate Eq. ( 2.40| ) into diagonal and off-diagonal parts, using the Hermiticity of the density matrix p: 

k>l k>l 

1 — ^ ^ ^mm^kkPkk "t - ^ ^ (Smm,kl ^mm^lk )Re(p k i) + i ^ ( S "unM - S mm jk)lm(p k i). (2.42) 

m,k ra,k,l m,k,l 

The first term of Eq. ( [2.42 ) is always 1 because of Eq. ( |2.41 ). If we require Eq. (2.42) when the off-diagonal terms 
in p arc p kl = S k ,k Si,i (k > I), we obtain 

y ' (^mm.tolo +<5 mm ./ fc ) = 0, (2.43) 
m 

and setting the off-diagonal terms in p to p k i = i5k,k $i,i {k > I) gives 

*^^(Smm,k la ~ S mm .l k ) = 0, (2.44) 
m 

Adding these equations, we obtain the desired result 

y^£mm,fc i =0, k Q ^ l . (2-45) 

m 

□ 



D. Perturbation theory 



In this section we develop a perturbation thoery in the coupling A for the superoperator S\j- The calculation 
will assume the diagonalizability of S\ t t- If all the eigenvalues of a matrix M are distinct, M is diagonalizable 
Therefore in many cases of interest for equilibration, this assumption for S\j will be correct. An example of a simple 
superoperator that is nondiagonalizable is the following. The superoperator S operates on BiTi^) and is given by 
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s(N)O'l) = o, 

5(|1)(1|) = |2)(2|, 

5(|2)(2|) = |2)(2|, ^ 
5(|3)(3|) = |1)(1|. 

The eigenvectors of iS are for all i ^ j, the state |2)(2| and |1)(1| — |2)(2|. This example shows that nondiago- 
nalizability is not a property particular to superoperators describing quantum operations but is also found in classical 
Markov processes. 

One can formally expand the superoperator Sx,t as a power series in the coupling parameter A, 

Sx, t = 5 t (0) + AS t (1) + A 2 sf ' + A 3 5 t (3) + . . . . (2.47) 

In section [IE we will explicitly calculate the expressions for these expansion operators. We will show (Eqs. ( |2.70| )- 
( 2.73 )) that condition Eq. (2J)) implies that is zero for all t. On the basis of this expansion, we will make a 
perturbative expansion of the eigenvalues and eigenvectors of 

M = ^°)+A M W+AV 2) + --., (2.48) 
X = X (0) + Ax (1) + A 2 X (2) + . . . . (2.49) 

Assuming that the perturbation expansion exists for this non-Hermitian operator, it will have the same structure as 
in the well established procedures familiar in quantum theory for bounded Hermitian operators (see textbooks on 
quantum mechanics such as J30| ] or J3l[ | for a more mathematical background). 
In the representation of Eq. (2.24) reads 

(sl 0} ) mnM = (t/'Wt/'V, (2.50) 

where U = e lHs . Unitarity of , as a matrix operator on , follows from 

'^2(St° ) )mn,kl(St° ) )kl,ij = ^ (U* )mk (U* ^ )/„ (U* )jl (U* ^ )ki = S m iSj n . (2-51) 

The eigenbasis of 5 t is formed by the set of matrices p nm = \n)(m\ where \n) are the eigenvectors of H s . These 
eigenvectors come with eigenvalues ^i°n m ' 

r„ „(0) _ At(B„-E m )xN,N (2 r 2] 

where E n are the eigenvalues of H s . Thus all density matrices of the form p nn , and mixtures of these, have degenerate 
eigenvalues /4°2n = 1- ^ the spectrum of H s is non-degenerate (we assumed this in section [II B| ), then all other 
eigenvectors p nm for n ^ m have non-degenerate eigenvalues. These eigenvectors p nm form an orthonormal set with 
the vector inner product on C N , 

Tr PnmPkl = dnkSrnl- (2-53) 

To carry out the perturbation theory, we switch to a ket notation for the density operators and a matrix notation for 
the superoperators. This will make it easier for us to perform the necessary manipulations of degenerate perturbation 
theory, in which the degenerate sector is isolated and a diagonalization performed within it. 

We first organize the diagonal, degenerate part of this vector space to be indexed. To be specific, we introduce an 
orthogonal basis in this vector space such that 

\<pf ] )=Pu, 1<%<N, (2.54) 
\^f m ,n))=Pmn, 1 < m, n < N, m + n. (2.55) 

In the second equation the indexing i can be made consecutive by choosing 

i(m, n) = nN + m — ^n(n + 1), m > n, , , 

i(m, n) = ^N(N — 1) + mN + n — ^m(m + 1), n > m. 
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This organizes this new vector space into a direct-sum form C N = Cd®Cnd, where 11 D" and "NB" stand for diagonal 
and nondiagonal (or, degenerate and nondegenerate). Cd has dimension N and Cnd has dimension TV 2 — N. 

From the discussion above, we note that the degeneracy is lifted in lowest order by the second-order part of the 

(2) (2) 

superoperator S in the D sector, which we will denote <Sp D . Assume that D is diagonalizable via the similarity 
transformation 

MS^oM- 1 (2.57) 

where Sjj D is a diagonal matrix (the tilde will denote quantities expressed in the new basis Mjj © 1jvd|0^)j which 
is in general non-orthogonal). In this new basis, the degeneracy of the diagonal terms of S is lifted to second order in 

A (the diagonal terms can be written to second order as /x« = 1 + X 2 S^'), and since the largest off-diagonal terms in 
the D sector are now third order, given by 

X 3 MS^ ) D M- 1 = X 3 S^ D , (2.58) 

the condition for the successful application of non-degenerate perturbation theory is now satisfied, assuming that no 
additional, accidental degeneracy occurs. (The condition is satisfied from the start in the ND sector.) Its form is 
essentially no different from the conventional perturbation expansion | f30| for Hermitian operators. This expansion 
for the eigenvalues is 

Mi = M f+A 2 «S^ + 0(A 3 ). (2.59) 

The form of this expansion is different depending on whether i 6 D or i € ND, but only at 0(A 4 ). The perturbation 
expansions for the eigenvectors are 

_ 5 (3) 

jeD, j^i s u s jj 
5 (2) 

|^) = l0f)+A 2 El0f) (0) r ' (0 ) + Q ( a3 )' ^ ND - ( 2 - 61 ) 
This expansion indicates that there is no mixing between the D and ND sector s until second order in A. This expansion 



strategy will be taken up again in the numerical simulations, Sec. II G (Eq. ( 2.116 )) 



We note that thi s sep aration of the superoperator into D and ND sectors permits us to write the action of the 



superoperator Eq. ( |2.47 ) using a more informative expression in which these sectors are almost decoupled: 

(S\,t{p))nn = ^ Pnm,t Pmrn + $Pnn, (2.62) 
m 

P'nm.t &nm © A («Sj *)nn,mrn ~t~ A (<S^ ^)nn,mm "t" • ■ ■ , (2.63) 

Pmn=\ 2 (St 2) )nnMPkl+\ 3 ^ (^UnMPkl + ■ ■ ■ ■ (2-64) 

k,l,k^l k,l,kjtl 

Note from Proposition [|that P nm ,t is exactly a stochastic matrix; therefore the dynamics in the D s ector is that of a 



classical Markov process, up to second order in A (since the contribution from the ND sector, Eq. ( 2.64 ), is 0(A 2 )). 
The dynamics inside the ND sector is also simple: 

(Sx,t(p))nm = Mt°nm Pnm + A 2 ^ ( S t )nm,klPkl + ■ ■ • , n ^ TO. (2.65) 

k,l,k^l 

So, to (9(A 2 ), there are no contributions to this equation from the D sector. So, the low-order dynamics in the ND 
sector simply involves a scala r mu ltiplic ation of the off-diagonal components of the input matrix p. 

The simplications of Eqs. ( [2.62 ) and ( 2.65| ) makes it possible to answer questions about the uniqueness of the fixed 



point and, in principle, the mixing properties of a repeated application of S\,t : using techniques from classical Markov 
processes [B2|. The splitting in two sectors, each having its own relaxation times, is similar to the phenomenological 
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description of a relaxation process by means of Bloch equations or the Redfield equation 1 33 1 . This description in 
terms of the longitudinal relaxation time T± (D sector) and transversal relaxation time Ti (ND sector) is, for example, 
used in NMR |33) . 

Of course, the "smallness" of the operators X 2 S^ , A 3 5^ 3 ^ , . . . compared to will determine how fast the pertur- 
bation series converges. We will calculate the eigenvectors of <Sa,* to zeroth order in A and the eigenvalues to second 
order in A. The stochastic matrix P nm ,t is determined in this approximation. The justification of this approximation 
will be given when we explicitly determine the expressions for in section HE, where we set bounds on A and t 
such that indeed A 2 and higher order corrections are small within some norm (for example the | . | o given m ||o4|,p5|| ) . 



E. Calculation of expressions 



Here we will calculate the elements of the sup erope rator described in the last section to lowest non-trivial order in 
A (A 2 ). Truncating the expression for P in Eq. (2.63) to second order, Q nm ,t is defined by the expression 



Pnm.t ~ $nm ~t~ A Qnm.t- 



And taking p in Eqs. ( [2.65 ,2.64) to second order, and using Eq. (2.52), we define v n m,t by 

Vn m ,t~e lt{E ~- E ™\l + \ 2 v nm , t ). 



(2.66) 



(2.67) 



In this section we will find expressions for Q nm ,t and v n m,t and exhibit the regime in which they give a valid description 
of S\ : t- We also show that for a large enough bath, the equilibrium state is the fixed point of the map S\,t- We 
discuss under what conditions this fixed point is unique. 

We will use operators in the Heisenberg representation. We denote such operators (for example on the system) as 



A t = e lH ^Ae 



iH.t 



The total Liouvillian C is defined as 



-iCt 



{p®Pb,p) = U t (p®p b ,p)U' 



tt 



(2.68) 



(2.69) 



One can expand the operator e in a perturbation series in A |2C(], take a partial trace over the bath and identify 



the operators 5 t (0) = e~ iCst , S t (1) and S t (2) in Eq. (|2~47|) 



(i) 



iTi b / dt'e- i ^ +c ^ t - t '\C sb e- i ^ +c ^ t ' , 



and 



S, 



(2) 



s, 



Tr, tdt' f dt" e -^+^X*-*')£ -i(£.+£.)(t'-t") £ -<(£.+£,)*" 



(2.70) 



(2.71) 



(2.72) 



where pt> is the time-evolved (with H s ) p and p b ,/3 t , is the time-evolved (with H b ) p b ,p. The equilibrium state p bi p is 
invariant under unitary evolution with e tHbt and thus pb,/3 , = Pb,0- We then use the cyclic permutation invariance 
of the trace and H s b = S ® B to rewrite equation (2.72) as a simpler sum of two terms 



First we consider <sf . We use Eq. ( |2.69|) and Eq. (|2.14| ) to rewrite S\ L> acting on pig) p b ^ as 



S^ip® Pbi0 ) = -i\Tr b / dfe iff .C*-0 ® e ^(*-0 [H sb ,p t , ® p b a, ] e~ iH ° ('"'') ® e^C*"*'), 



S^ip® Pbs) = -<A / dt' e iH ^-^S Pv e- iH ^-^ - e iH ° pfSe~ iH ° <*"*') 
Jo L 

Then the condition Eq. (2.9) implies that Sj^fa® p bi p) is for any p. 



Tr b Bp bt 



(2.73) 



Let us consider the second order term. The expression for S} reads 
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s (2) = _ e ~iC s t 



f dtf [ dt" {h{t' -t")S-t>S- t "p-h{t" -t')S-fpS- 
Jo Jo 

h(t' - t")S-t»pS-t> + h(t" - t') P S-t>>S- t i) , 



where h(t) is defined as (BBf)^. We write 



h(t) = / duo e ltu) h(u>). 



(2.74) 



(2.75) 



Let S nm be the matrix elements of the interaction S in this eigenbasis of H s , S nm = (n\S\m). Now we can find the 
express ion fo r Q m n.t — {S\ 2 ^) m m,nn- From Eq. ( |2.74|) after integration over the variables t' and t" and with the use 
of Eq. j2~75|), we find: 



Q i, 



du> h(oj) 



\S mn \ 2 {l - coat{uj - E n + E m )) v-^ 6 nm \S nl \ 2 (l - cost(w - E n + £?,)) 



(lu - E n + E m ) 2 

For the "decay factor" v nm ,t in the ND sector we find 

/*°° j It ,\^S nn S mm {l-costw) 
v nm ,t = do; h(u!) 



{W-En + Eif 



f(t,u,E n )- f*(t,uj,E m ) 



with /* the complex conjugate of /. The function / is given by 

\Si n \ 2 (l-cost{uj-E n + E l )) 



Bef(t,u,E n ) = ^2- 



and 



Imf(t,uj,E n ) = J2 



\Sif, 



{oj-En+Etf 

svcitloj — E n + Ei) 



uj- E n + Ei 



1 - 



t{u-E n +Ei) 



(2.76) 



(2.77) 



(2.78) 



(2.79) 



We will now look at the idealized case, i.e., we take the limits (remember k is the number of qubits in the bath) 

(2.80) 



P nm ,\ 2 t = km lim P nm ,t, H r , 

t— *oo,A— 1-0 k — >QO 
constant A^t 



X 2 t = e^-E™) lim l im (i + \ 2 Vnm +). 

t— »oo;,A— *0 k — »OG 
constant t 



When the bath is infinitely large, it will have a continuous spectrum; h(u>) will be a smooth function. The rate of 
interaction vanishes, but as we take the limit t — > oo, there is an effective non-zero interaction that is proportional to 
\ 2 t. Recall that 



8 (x) = lim 



1 — cos(te) 



t— «x> tnx 2 



(2.81) 



where S(x) is the Dirac delta function, which is defined as f_ dxS(x) = 1 and, Vi ^ 0, S(x) = 0. With the use of 
the S function we find 



Pmn,XH = <W1 - X 2 t2nJ2\Sni\ 2 h(En ~ E t )) + \ 2 t2v\S mn \ 2 h(E n - E m ), 



and 



with 



P J nm,\ 2 t 



= e^ E "- E ^ (l + X 2 t2TrS nn S mm h(0) - \ 2 tng(E n ) - \ 2 tng*{E m ) 
Re g(E n ) = ^ \S ln \ 2 h{E n - £?,), 



(2.82) 



(2.83) 



(2.84) 



and 
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Im g(E n ) = V ^ duh(w) £ ^ - (2.85) 
where "P is the principal value of the integral. In order to see in what regime the perturbation theory is correct, 



we check whether the process descr ibed by Eq. (2.82) and Eq. (2.83) corresponds to that of a TCP map. First 
we verify Property in Eq. (2.83); the eigenvalues of |n)(m| and |m)(n| are related by complex conjugation, or 



tKim x 2 t ~ t 1 mn.x 2 t- The trace-preserving property (also in |lj) is also obeyed: 

Y,P mn ,xn = l- (2.86) 

m 

Complete positivity of the map implies that P rnn ,x 2 t must be a matrix of probabilities, that is we must have P m n,\ 2 t > 
0. Thus the first necessary condition for the validity of the perturbative approximation is 



Condition 1: Vn: A 2 K 



Eq. ( 2.86| ) and Eq. ( 2.87 ) together ensure that P m n,x 2 t is a stochastic matrix. Complete positivity also implies via 



Proposition ^ that |/i nmj A 2 tl !• in order that |1 + X 2 ta\ < 1, where a is some complex number, we must have that 



Re a < and X 2 t < 2/|Re a\. This real part in Eq.(2.83) is indeed negative as h(ui) is positive, and we obtain a new 
condition: 

Condition 2: V m, n : X 2 t « -^-—^ ——^J^———^ ^ ]Sim?h(Em _ Ei)] ■ (2.88) 



Note that this condition is quite similar to the condition in Eq. ( 2.87 ) 



It is not hard to see that the stochastic matrix P m n,x 2 t obeys detailed balance for the equilibrium distribution: 

P m n,XHer pEn = P nm .x 2 te- pEm , (2.89) 

as the equilibrium condition of the bath implies that 

Ii(-lu) = e- 0u h(oj). (2.90) 

Thus the equilibrium density matrix p s ^ is a fixed point of the idealized equilibration process. To consider whether 
this fixed point is unique, we note the following: If a stochastic matrix M is such that all its matrix elements My > 0, 
then M has a unique eigenvalue equal to 1 ffij. If Condition 1 is obeyed, we indeed have P m n,x 2 t > and therefore the 
absolute value of the second largest eigenvalue (in the diagonal sector) is smaller than 1. For the off-diagonal sector, 
Condition 2 says that the largest eigenvalue in the off-diagonal sector is strictly smaller than 1 in absolute value. Thus 
under these conditions, with Proposition we can conclude that the process converges to the equilibrium state. The 
expression of P mn ,x 2 t coincides with the derivation given by Davies |25| for small \ 2 t. 

One can help to speed up the process in the off-diagonal sector by "dephasing" ; that is, after having the system 
and the bath interact for some time t, we perform the operation 

V a (Ps) = -Ye m ° s p s e- lH ° s , (2.91) 



a 

s=Q 



which can be implemented with the assistance of an extra register in the state ^= X)s=o I s ) which is used to condition 
the evolution U = e lHsS and subsequently traced out. The dephasing has the effect of canceling off-diagonal terms in 
the eigenbasis of the system, i.e. 

lm^I? Q ^ QiklPkl =^2a kk p kk . (2.92) 
\ fcJ / k 



A complete dephasing can in general not be achieved in polynomial time in n (see section III), and thus must be 
understood as an extra aid but not a solution to the equilibration problem. 

From the expressions for P mn ,x 2 t and Pnm,x 2 t we can understand the physical picture of the interaction between 
bath and system. The system makes a transition from (eigen) level n to level m, when (1) S mn is non-zero, (2) the 
bath is capable of "receiving" this quantum of energy AE = \E m — E n \, that is, it has a matching energy difference 
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\u>i — — AE and (3) By is non-zero. Furthermore, the more such transitions there are, the faster the off-diagonal 
matrix elements decay. This confirms the intuitive picture that one might have of equilibration. Note also the 
similarity with the Fermi Golden Rule p2|]3l| that describes the transition probability from eigenlevel n to m in a 
unitary evolution that is perturbed by a time-dependent Hamiltonian. 
For a finite-dimensional bath, we can express h(t) = (BB t ) as 

h(t) =J2^ u "'-^\B kl \ 2 e-^ k /Z b , (2.93) 
fc,; 

where Bki = (ki,\B\lb) with being the eigenstates of the bath Hamiltonian Hi, and Zb the partition function of the 
bath. Taking the limits t — > oo and A — > before letting the bath grow large leads to divergent expressions for P mn ^ t 
and fJ>nm,\ 2 ti suggesting that the perturbation theory fails in this regime. This is not surprising, as the fmiteness of 
the bath together with the limit t — ► oo will lead to Poincare recurrences (only the interaction cycle time is long due 
to A -> 0). 



F. The inverse quantum Zeno effect 



In our numerical studies (sections II G and II H ) we have observed a phenomenon that one might call the inverse 
quantum Zeno effect. It is a way of mapping an arbitrary initial state onto the completely mixed state Iff by 
interacting repeatedly and strongly with the state for a very short time. Here we will give a theoretical analysis that 
expl ains t his observation. Consider the weak coupling expansion Sx,t 
Eq. ( 2.71 ). We expand these operators around t = 0: 



5 t (0) + A 2 ^ 2) + C(A 3 ) with S ( t z> given as in 



:(2) 



(2) 



2\2 



t 2 X 



SxAp) =P~ U[H S , p] + ~-£-(\Sp, S] + [S, P S}) (B 2 ) b + 0(t 2 , X 3 t 3 ). 



(2.94) 



In the limit A — > oo, but t — > 0, and constant A t, the higher order terms 0{t , A t 3 ) will vanish. Thus we see that 
the fixed point of in this limit (assuming non-zero (B 2 }b) must obey 



[H s ,p} = 0k[[S,p],S}=0. 



(2.95) 

Notice that if we take the differential form of Eq. ( [j.94| ) and the prescribed limit, the equation is of the Lindblad 
form, Eq. ( 2.17 ). The state In certainly meets the requirements of Eq. ( 2.95| ), but is it unique? If S and H s are such 
that they have no eigenspaces (except for the full spac e) in common, and both have a non-degenerate spectrum, we 
can show that ljy is the unique eigenvector. Eq. ( 2.95 ) requires that either [S, p] = or [S, p] is diagonal in the same 
basis as S. If [S, p] = but also [H s ,p] = 0, then p can only be the state ljv. What happens if [S, p] is just diagonal 
in the same basis as S7 Let \n) be an eigenvector of S with eigenvalue X n . We have for n^ra 



Rewriting this expression gives 



(n\[S,p]\m) = 0. 



\/n,m,n^m (n\ p |m)(A n — A m ) = 0. 



(2.96) 



(2.97) 



Now, because p is diagonal in the basis of H s as [H s , p] — and H s and S have no eigenvectors in common, there 
exist n and m such that (n\p\m) ^ 0. But the eigenvalues of S were non-degenerate, thus we obtain a contradiction. 
□ 

When ljv is the unique eigenvector of this process, then, with the use of Proposition [| the repeated application as 
in step 4 of the Equilibration algorithm I will eventually bring the system to the state ljy. 

We showed that for this "inverse quantum Zeno" effect to occur S and H s have to be such that they have no partial 
eigenspace in common and both have a non-degenerate spectrum. If we assume that S and H s are c-local with c 
larger or equal to 4, then this does not impose a very strong constraint on S and H s : the effect will occur for a generic 
S and H s . 
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G. Specifications of the numerical simulation 



The main purpose of this study is to understand the effects of bath size and the choice of bath and interaction 
Hamiltonians for a specific system Hamiltonian. In Table | we list some of the choices that have been made in the 
num erica l ana lys is. We have randomly generated the elementary Hamiltonians hi that make up H s ,Hb and H s b, 



Eqs.(2.5), (2.6), (2.7), with a measure A4. We choose the diagonal elements of each hi uniformly in [—a, a], where a 
is sampling scale in Table Q. The absolute value of the above-the-diagonal elements of a hi are chosen uniformly in 
[0, a] and its phase is chosen uniformly in [0, 2ir]. The below-the-diagonal elements of hi follow from Hcrmiticity. This 
defines M.. Note that M. is not a unitarily invariant measure. 

We take the Hamiltonians S and B as sums of all possible local 2-qubit interactions (c s = 4 in Table |) . For the 
Hamiltonian of the system H s we also take a sum of all possible local 2-qubit interactions. Note that this includes a 
set of Hamiltonians that exhibit frustration, for which we don't expect equilibration to be particularly fast. 



In section II E we observed that matching energy differences between bath and system are an important ingredient 
in the equili bration of the system, which is consistent with the intuitive picture of equilibration that was sketched in 
section II A. However, as we do not know the eigenvalues of the system, we can only pick our bath so as to optimize 
the chance for matching level differences. The sampling scale of the bath f(n,k,c s ,Cb) is determined by roughly 
optimizing these coincidences, AEf, = AE S . 

Consider the density of states p s (E,a s ) of the system (the distribution of eigenvalues generated by the measure 
M.) and the density of states Pb(E, a,b) of the bath. Here a s is the sampling scale of the system which we set to 1 (see 

Table |). The quantity [Tri? s ]>i is the mean and ^- TlH ^ M is the variance of the distribution p s (E , a s ). The choice for 
M. ensures that the distributions are symmetric around E = 0: 

[Tr H S ] M = [Tr H b ] M = 0. (2.98) 

To optimize for matching we choose the variances to be equal: 

{ b\M 



[Tr H*] M _ [Tr H$ 
N K 



(2.99) 



For large K the bath distribution will be Gaussian (central limit theorem), whereas the system distribution will be 
similar to a Gaussian distribution for large N (see Fig. ||) . Thus, setting the variances equal brings the distributions 
close together. 

Consider first [Tr H 2 ]m- It is straightforward to calculate the variance of the eigenvalues of a qubit bath. Given a 
2x2 Hermitian matrix m^-, the eigenvalues e± = \{pii\ + ^22 ± yj (m,\\ — TO22) 2 + 4|mi2| 2 ) have 

1 f ab f ab f ab la? 

\ e V\M = -r^ dm n dm 22 d\m 12 \ e\ = — (2.100) 
4a b J-a h J~a b Jo 

Let Vi be some ± pattern i of length k, corresponding to selecting e+ or e_ for each qubit bath. Let E Vi be an 
eigenvalue of the full bath, i.e., E Vi = Y^m=i e vi[m] where Vi[m] indicates that we select the mth bit in Uj. Then 

i—l 



We calculate [Tr H^] M = J2% j[\{H s )ij\ 2 ]M f° r n > 2. We can write 

(?) 

£[|(ff s )tfl 2 U =E EtKMofU (2-102) 

i,j i : j m—1 

where h m is the mth local interaction Hamiltonian. We have used [(h^.)ij (h m )ij]M = 0- Each row of h m has only four 
non-zero entries as the dimension of the local Hamiltonians c s was set to four. Using the fact that [\{h m )ij \ 2 ]m = 3 
for all interaction terms m, we obtain 



[Tr H 2 S ] M 4 



n 



N 3 V 2 

For n = 1, we have E-Ej^tl^ = |, Comparing Eqs. (2.101) and ( [2.103 ) gives the expression for ab- 



(2.103) 
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a b = jf(n,M,2) 



(2.104) 



'l/k. Fig. |l| illustrates how this setting determines the density of states of bath and 



For n = 1, /(l,fc,4,2) = 
system. 

The numerical work consists of a calculation of the fixed point of 5a, t as a function of t for a fixed A and the second 
largest eigenvalue for different b aths and different systems and temperatures. We follow a numerical procedure based 
on perturbation theory (Section |IID| ) to perform a stable numerical evaluation of these quantities. We can trust the 
answers from the numerical procedure only if we are in the regime in which perturbation theory is correct. This 
regime was heralded by the two conditions Eq. (2.87) and Eq. (2.88) in section O. Whether these conditions are 
obeyed depends on the specific choices of H s , Hb and S and B. We prefer to reformulate these conditions here such 
that they are obeyed for the average bath, system and interaction Hamiltonian o btain ed by sampling using M. and 
the sampling scale. As the conditions are very similar, we take the first one, Eq. (2.87), and reformulate it as 



c(t) = \ 2 t2ir 



NK[S 2 ] m [B 2 ]m 



W b 



< 1. 



where [S* 2 ]^, the average matrix element, is defined as 



^ = ^ £[l<s«l 2 u = wl^u, 



and similarly for [B 2 ]m- Wb is the spectral width of the bath, i.e., 



Wd = 



[Tr H 2 ] M 
K 



Here we indicate the approximations made in obtaining Eq. (2.105) from Condition 1 (Eq. (2.87)) 

\H2-K^\S ln \ 2 h{E n -E l )^l. 



(2.105) 



(2.106) 



(2.107) 



(2.108) 



Using Eq. ( 2.93| ) and Eq. ( 2.75| ) we write the h function as 

h(E n -E l )=Y d S((E n - Ei) - (uj k - uj m ))\B km \ 2 e-^/Z. 



k.m 



(2.109) 



We will approximate the matrix elements |-Bfc;| 2 as constants and replace them by their average [£? 2 ]x. Then we can 
use density-of-states arguments to approximate the m sum over the S functions by the inverse of the average spacing 
between the 5 functions; this spacing is given by W b /K: 



S((E n - El) - (0J k - U} m )) 



K 
Wb" 



(2.110) 



With these approximations, the partition-function sum over k in Eq. (2.109) becomes exactly one. So, Eq. ( 2.109] ) 
becomes 



h{E n -Ei) 



K[B 2 \ M 
W b 



Now Eq. (2.108) is 



X 2 t 2ir 



K[B 2 ] M 
W b 



< i. 



(2.111) 



(2.112) 



If we again approxima te the matrix elements |Sz n | 2 as cons tants and replace them by their average [S^Jx, an d note 
that the I sum in Eq. (2.112) has N terms, we obtain Eq. ( 2.105| ). 

For the simulations we have performed, we can find the values for [S 2 ]^ and [£? 2 ].m (note that these Hamiltonians 
have locality parameter c — 4, as does the system Hamiltonian H s ) and obtain the expression 
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< 1. 



for n > 1 and k > 1. For a qubit system, n = 1, and k > 1 we obtain 



ci(t) = \H 



_ - 8ttV2 



3%/3 



< 1. 



(2.113) 



(2.114) 



The quantity c(t) in Eq. ( |2.105 ) will function as a rescaled time which depends on the strength of A and the size of 
system and bath. In the regime where c(t) < 1 we expect a pertubative calculation of the eigenvectors and eigenvalues 
of the superoperator to be fairly accurate. The dimensionless parameter associated with the temperature is given by 



I3' = (3W S 



(2.115) 



where W s is the spectral width of the system, Eq. (2.107) (W s — Wb). From here on, (3 w ill re fer to this scaled 
dimensionless parameter. Instead of expanding the superoperator S in a series in A as in Eq. (2.47), we write 



s. 



s, 



(0) 



(2.116) 



~(2) 

where all higher order terms are grouped in S t . The calculation of eigenvalues and eigenvectors then follows the 
analysis of Section II D. We find that the choice for the bath and the interaction Hamiltonian influences whether the 
equilibration will succeed or not. Let 



V 



Ps,/3 - P0 \\tr, 



(2.117) 



where po is the unit eigenvector obtained from the numerics. In Figs. ^| and ^ two extrema in dynamics are shown, 
each corresponding to a different choice for the system, bath, and interaction. In Fig. |^ the equilibration is successful, 
whereas in Fig. [| the equilibration fails. Rd is defined as 



Rn = 



1-\k d \ 

m 



(2.118) 



where kd is the second largest eigenvalue in the diagonal sector and c(t) is the average coupling strength in the time 
interval that we consider, which is c(t) € [0, 0.3] here. Similarly, we define 



Rnd 



\knd\ 



(2.119) 



for the nondiagonal sector. 



H. Numerical results for equilibration 



We are interested in how well a randomly chosen bath and interaction equilibrate a system and how these averages 
are improved by choosing larger baths. As the mixing rates and the distance to the equilibrium state will in general 
be oscillating functions of the scaled time c(t) (see Fig. ||) we will compute time averaged rates over a reasonable 
interval in c(t), 

[c(t init ) =0,c(t end ) =0.5], (2.120) 



such that we are in the realm where perturbation theory is valid, Eq. ( 2.105| ). We denoted these time averages (not 
to be confused with bath averages) as Rd and V for the time averaged trace distance, Eq. (2.117), etc. In Fig. [| we 
present histograms that show how, for a given fixed system and interaction, the equilibration process is different for a 
set of randomly chosen baths with fixed dimension. The insets show the distribution for the lowest bin. The vertical 
axis denotes the percentage of baths (the interval [0%, 100%] is given as the interval [0, 1]) for a certain distance and 
rate. We observe that the diagonal rate distribution is very broad, and therefore the mean of the distribution is not 
a very good (or a very stable) measure of the generic behavior. Furthermore, we find that the rate in the diagonal 
sector is much worse than in the nondiagonal sector and thus is the dominant factor in setting the mixing time. This 
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conforms to the pattern in many quantum systems, for example for nuclear spins as observed by NMR, for which T\ 
is generically larger than T2 33 1. 

To study the dependence on f3 and on the dimension of the bath versus the dimension of the system, we compute 
the following data. We pick a system Hamiltonian H s of n qubits that has some well spread out spectrum. We set 
the dimension of the bath and then we randomly pick both the bath Hamiltonian and the interaction Hamiltonian. 
Means are denoted as [.]Mb- For the rates we look both at the mean and the median. The median is denoted as 
[[.]]jki, see Fig. ||. The results for n = 1, 2, 3 and 4 are shown in Figs. We have given the median when the mean 
does not give a good representation of the distribution. 

These data clearly indicate that larger baths improve the process of equilibration, both in the rates (D and ND) 
as well as in the closeness to the equilibrium state. The effects are the most pronounced at low temperature, where 
equilibration is in general harder as the system must relax to a single pure ground state. To understand the closeness 
scale, we show in Appendix [A] how far apart two arbitrarily chosen density matrices are; this number lies around 1 
for the dimensions that were considered. For these estimates, we see a trend towards approximations getting worse 
for larger system sizes for low temperature. The scaled rates [RD\M b and [Rnd]mi, seem to be fairly constant, thus 
we see behavior that suggests that the rates are polynomially related to both system and bath number of qubits. We 
also observe that the nondiagonal rate (ND) is always higher than the diagonal rate (D). The data show a system 
Hamiltonian dependence, that is, the average equilibration for n = 4 seems to be more succesful than for n = 3. We 
also observe that the difference between T\ and Ti becomes smaller with increasing (3 (lower temperature). Thus, in 
conclusion, it seems if we pick a bath size (in number of qubits) that is polynomially related to the system size (note 
that the number of eigenvalues is then exponentially related), the rates of relaxation are polynomially related to the 
system size (in qubits); however the relaxed state could be still fairly far away from the true equilibrium state for 
large system sizes. 



III. EQUILIBRATION II 

We present an alternative to the algorithm in section H. This algorithm relies on the technique for the estimation 
of eigenvalues, originally given in Q] (see J37||3l|l). This eigenvalue estimation routine has also been used as a building 
block in an interesting quantum algorithm in |p_9|| and p9| . 

Let H s be the c-local Hamiltonian with non-degenerate eigenvalues as in section [n| Order the eigenvalues as 
Eq > Ei > . . . > En. 

Definition 2 Equilibration algorithm II. 

1. Initialize the system in the (infinite temperature) completely mixed state In- Also add one m-qubit register set 
to |00...00)(00...00|. 

2. Compute eigenvalues with the use of the Fourier transform and dephase in computational eigenvalue basis, 
which will result in state 

2V-1 2 m -l 

J2 Yl P(s,n)\n)(n\(g>\s)(s\, (3.1) 

n=0 s-1 

where p(s, n) is a distribution, peaked at s ~ E n for large m. The dephasing is a simple superoperator T> on the 
eigenvalue register that operates as 

V{\si){ Si \) = \ Si )(si\, V(\ Si )( Sj \) = 0. (3.2) 

3. Prepare an additional N -dimensional quantum system, the bath, also in ljv- Add a m-qubit register and one 
qubit register set to |00 . . . 00) (00 . . . 00|. 

4- Compute eigenvalues of the bath as for the system in step 2. 

5. Interact system and bath according the following rule TZ ("partial swap"): 

( \m,n)\s,t)\0) ift<s 
Un\n,m)\a,t)\Q) = < p/2 f J (3.3) 

[iPst |M)|0) + yl -p p st \s,t)\l))\s,t) ift>s 

where p^ t = e~^^~ s \ 
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6. Trace over the single-qubit register, all bath registers, and the eigenvalue register of the system. The system 
will be in some state 



Ps 



E 



\n)(n\ 



(3.4) 



The steps 2-6 define a TCP map S, iS(ljv) = p s . 
7. Repeat steps 2-6 r times such that 

|| 5 r+1 (|000 . . . 00) (000 . . . 00 1) - 5 r (|000 . . . 00) (000 . . . 00 1) || tr < e, 



(3.5) 



for all r > ro and e is some accuracy. 

The advantage of this algorithm is its simplicity and its similarity to a classical algorithm; we create a Markov 
chain in the eigenbasis of the system. The disadavantage of the algorithm is that it is very likely to be slow; the 
computation of the eigenvalues to high accuracy with the use of the Fourier transfrom is very likely to be exponential 
in the number of qubits of the system and has to be performed twice, for system and bath, in each round of the chain. 
First, let us show that in the case when the eigenvalues are computed exactly in steps 2 and 4, i.e, p(s, n) = SE n , s /N 
the Markov chain equilibrates the system. Recall pTj that the routines of steps 2 and 4 compute rescaled eigenvalues 



E 'n — fiE n + fi 



(3.6) 



with /i and f% depending on the maximum and minimum eigenvalue (of which we assume that we can find an estimate) 
such that E' n £ [0,27r). In the following we will drop these primes. The chain that is created can be represented as 



,(*) 



n)(n| 



(3.7) 



where a 



(*) 



(fc-i) 



Note that £ ro P n 



P n ^m- We have 



Pn- 



l_ 

N a 

&(i + £*<«(i-j&)) i£E m = E, 

J_„/3 

N^nm 



if E m < E n 

if E m = E n 
if E m > E n 



1 as required. The equilibrium state Eq. (1.1) obeys the detailed balance condition: 



Vn, m P n 



-0E n _ p 



-@Er, 



(3.8) 



(3.9) 



All the matrix elements of the Markov matrix P n ^ m are nonzero. Therefore the chain will have a unique fixed 
point which is equal to the equilibrium state due to detailed balance. Thus for all probability distributions a n we 
have 



lim V a„pW TO 



-I3E„ 



(3.10) 



Notice that it is not hard to prepare the initial states of system and bath. One way to make the completely mixed 
state ljv is to make a maximally entangled state ^= ^2^=^ an d trace over the second register. This takes 0(n) 
steps. The partial swap in step 5 can be implemented with 0(n) elementary qubit steps. The dephasing in step 2 is 
introduced to keep the form of the algorithm clean, but it does not affect its output. This dephasing is implemented 
by measuring the eigenvalue register in the computational basis and discarding its answer. When using an m-bit 
eigenvalue register the joint probability p(n, s) in the first round (after step 2) is equal to 



P(n,s) = — 



e il(E n -2Trs/2 m ) 

1=0 



(3.11) 



When p(n, s) is not a delta function on the eigenvalue, the Markov chain will still be in the eigenbasis of the system; 
It will be a concatenation of chains; the transition probability of this new chain is 
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J^Um = Z)K*|n)^-**P(»n|*). 



(3.12) 



where p(s\n) is a conditional probability, defined by p(n, s) = p(s\n)p(n), and P s ^t is the exact chain (when p(s\n) — 
Se„,s)- Note that J2 s P( s \ n ) = 1 an( ^ J2 m P( m \t) = ^> so tnat ^n->m ^ s a stochastic matrix. Let us make a few 
remarks about the behavior of such an approximate equilibration process. If this new Markov chain is close to the 
exact Markov chain, we can bound the deviation from the exact fixed point with perturbation theory pCfl . Let 



pi p 



(3.13) 



where E nm is a deviation matrix defined by Eq. ( |3.12 ). Let pa = p' s p ~ Ps,p where p' s g is the fixed point of the 
Markov chain PL m . Assume that P is diagonahzable. Let Y be the matrix defined as 



Y = (i _ P + pH)-i _ p(c 



(3.14) 



where p(°°) is the infinite iteration of P. Wc can write p(°°J = diag(l,0, . . . , 0) in the basis where the stationary 
state p Si( g is an eigenvector. In this basis, with diagonalizability, P is of the form diag(l, A2, . . . , Ajv)- We can then 
write 



Y = diag(0, 



1 



1 



l-K 1 



' 1 — A 



JY 



-), 



(3.15) 



where k is the second largest eigenvalue. For later use we note that the norm |||Y |||2 = h- K \ ■ It is possible to write 
the deviation pA in terms of Y and E: 



p A = (1 - YEY l YEp sJi 



(3.16) 



when E is small enough such that 1 — YE is invertible. This expression can be derived from P^'pa = 0, which 
follows from the uniqueness of the stationary state p. We now use 



PA \\tr< vN || p A H2, 



(3.17) 



as in Proposition |3j. Then using the expression for Y, Eq. ( 3.16 ) and Eq. ( 3.17 ) (see also below Eq. ( 2.37 )) we can 
bound 



PA \\tr< C n Tt p sj3 I 1 - — 



\E 



l-K V l-K 



\E 



l-K 



(3.18) 



and the rate of 



Thus the size of the correction pa will be determined by the strength of the perturbation \\\E 
convergence of the original Markov chain P. 

For a general H s , the computation of an m-bit approximation of the eigenvalues is likely to cost an exponential (in 
m) number of elementary gates. As there are 2™ eigenvalues, knowing the m — logpoly(n) bits of the values of E n 
still leaves groups of pjy^ eigenvalues indistinguishable. Thus only in very special cases, if the gates Uf can be 
implemented with a polynomial (in to) number of elementary steps (as in Shor's factoring algorithm H) is it possible 
to compute the eigenvalues to high accuracy efficiently. 

We have demonstrated a way to set up a Markov chain on a quantum computer that will converge to the equilibrium 
state for long enough time. For special Hamiltonians, there might be more efficient ways to tune and modify this 
kind of algorithm. The rule 1Z might be chosen to depend on other features of the eigenstates \n) and |m) as in the 
classical Metropolis algorithm where transitions are made between states that are related by local spin flips. There 
might be Hamiltonians for which the calculation of an eigenvalue, given the eigenvector, is efficient. Then there is 
the hard question of the (rapidly) mixing properties of the chain, that determines the computational efficiency of the 
algorithm. 



IV. (TIME-DEPENDENT) OBSERVABLES 

Given that we have prepared n qubits in the equilibrium state corresponding to a certain Hamiltonian H s , we 
can then proceed by experimenting and measuring. The simplest measurement that we could try to perform is the 
estimation of the expectation value of a c-local (Hermitian) observable O: 



21 



(O). = Tr p s<0 O. 



(4.1) 



As O is local, we write O — X)i=i Oi where each operator Oi operates on a Hilbert space of constant dimension c. 
We can calculate the eigenvectors and eigenvalues of each Oi rapidly on a (possibly) classical computer, which takes 
poly(rt, c) operations. If Oi has eigenvalues \Xi that are both smaller as well as larger than zero, we define Of as 

0+ = AO, + I mm Mfe |l) (4.2) 

max fc fi k + | mm fe p k \ k 

such that Of is positive semi-definite and has eigenvalues smaller than or equal than 1. If Oi has only positive 
or zero eigenvalues, we just "normalize" the operator by dividing by max^ p, k , and similarly if O has only negative 
eigenvalues. Let / be a positive operator valued measurement (POVM ^J) with operation elements A\^ and A 2 ^ 
and corresponding outcomes 1 and 2 such that 



4 

This measurement will give outcome 1 with probability 

p lti = Tr Of p etc. (4.4) 

The operators A\^ and A 2 ^ are given by 

Aij = C/ (diag +) 1 /2[/t an d^ = ;7 (l-diag +) 1 / 2 [/t, (4.5) 

where diag Q + is the diagonal form of Of and U a the diagonalizing matrix. We summarize these results in a Proposition: 

Proposition 5 The estimation of Tr pO where O is a c-local observable with precision 5 and error-probability e and 
p € B P os,i(Wat) (N — 2 n ) takes TO(ln i/<5 2 )poly(n, c) operations where T is the time to prepare the state p. 

Proof: All commuting observables Oi can be measured once for a sin gle preparation of p. To estimate a probability 
p with precision 5 and error probability e we need 0(ln \/8 2 ) samples f42f] . □. 

More interesting is an algorithm to estimate time-dependent expectation values. L et Q \ and Oi be two c-local 
observables. We consider how to estimate a time-dependent quantity (identical to Eq. ( l.lOj )) 



Tt p [O u O 2t ], (4.6) 

where Oi t is in the Heisenberg representation. Notice that C>2t, the time-evolved operator, will for general t not be 
local. Thus we cannot use Proposition ^. The way these quantities come about in linear response theory [Tl| ] provides 
the key for how to estimate them on a quantum computer. One considers a system that is perturbed at some initial 
time t = 0: its time evolution is generated by the perturbed Hamiltonian H s + \0\{t) {0\{t < 0) = 0) and A is small. 
After time t we consider the response of the system to the perturbation by measuring another observable O^. Notice 
that with Proposition |^, it is simple to perform this experiment. Linear response means that we take into account 
corrections of order A, but no higher order, in the estimation of 

6(0 2 ) s = Tr 2Pt - Tr O 2 p , (4.7) 

where pt is the time-evolved system density matrix. This first-order correction takes the form fl36| 

5(0 2 ) s ni\[ dt' Tr ^[Oi(*'),0 2i _ t ,]. (4.8) 



If the disturban ce \{t) — 0\8{t = 0) we find on the right hand side the correlation function of Eq. ( |4.6| ). The 
quantity of Eq. (|4.6|) is interesting, because it can be used to compute the simplest reponse of the system, the linear 
response of Eq. Q4.8| ), which we can directly estimate on our quantum computer, provided that both 0\ and 2 are 
local. But we are of course not restricted to a linear response regime: A is a parameter that we can tune freely. A 
sequence of measurements could determine higher response functions that will involve quantities such as 

(0 ltl 2t2 3ta ...O ktk ) s . (4.9) 
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V. CONCLUSION 



It seems that by asking the question of how fast real quantum systems equilibrate, we have opened a Pandora's box 
of hard-to-answer questions. If there are many simple quantum systems in nature that equilibrate slowly (that is, not 
in polynomial time) by any dynamics that does not require extensive preknowledge of the system, then it would be 
unreasonable to ask our quantum computer to perform this task efficiently. By relaxation in polynomial time we mean 
the following: in polynomial time we obtain a state that is within e trace distance of the equilibrium state where e is 
a small constant. It might be the case that leaving aside the classical phenomenon of frustration, relaxation does not 
take place in polynomial time. The idea here is that for a quantum system, the eigenbasis is not known beforehand, 
but must be singled out on the basis of an estimation of the eigenvalues, which is generically a hard problem. 

This however is not in contradiction with physical and experimental reality as we know it, as the quantities that are 
measured in an experimental setup usually involve operators on a small number of qubits; these are the experiments 
that can be done efficiently (in polynomial time) and thus do not necessarily probe the system's complete state. 
For example, the outcomes of the set of measurements tr^ <g> <Tj 2 ® . . . o~i n where <7j . is one of the Pauli matrices 
or 1, completely determines the state, but there are 4™ measurements in this set. In an experimental setup, we 
might randomly select a polynomial subset of them and there is some small chance of order pol ^ n ^ that these are the 
measurements that distinguish the equilibrated state from the present state in the lab that is supposed to approximate 
it. The estimates of time-dependent correlations could possibly be more sensitive to distance from equilibrium, as 
these involve time-evolved, non-local operations. The numerical study suggests that product baths whose size is 
polynomially related to the system can function as adequate baths in the sense of providing relaxation in polynomial 
time. The relaxed state could still be a rather rough approximation to the true equilibrium state, but, as we argued 
above, it might be a good starting point for subsequent measurements. 

We have taken the bath to be part of the (cost of) the quantum computer. In any experimental setup, there is a 
natural bath that is used to equilibrate and cool the quantum computer. Can we use this bath for a computational 
problem such as equilibration? Consider for example the NMR quantum computer ]l8| where computation takes 
place at room temperature. In the regime in which the heat bath has a non-Markovian character it has been shown 
to be possible to alter the Hamiltonian of the system and the coupling to the bath dynamically (see 0], but also 
standard books on NMR ||^| ) . These techniques could make it possible to simulate the time-evolution of a "designer" 
Hamiltonian and also to equilibrate the system to the equilibrium state of this designer Hamiltonian. 
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APPENDIX A: NORMS 



In this Appendix we give the definitions of several norms and inner products. The inner product between vectors 
in can be represented on B(Hn) as 

(X1IX2) = Tr x\x2- (Al) 

The trace norm l|34|.|35i is defined as 



|| A \\ tr = TrVlU. (A2) 
What makes this norm attractive is that it captures a measurable closeness of two density matrices p\ and p2 p5[ : 

|| Pi P2 ||tr= max]T \P?(]) - P A (j)\, (A3) 



A 

3 



where Pf 4 and P£ are the probability distributions over outcomes j that are obtained by measuring observable A on 
pi and pi- The matrix norm |||. |||2 is defined as 

||L4||| 2 = max ||Ar|| 2 . (A4) 
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where || . H2 is the Euclidean norm on C N : y/ffiv) for |t>) G C N . We have 

||Ar|| 2 <|||A||| a . Hzlh. (A5) 



In order to aid in the interpretation of the numerical results of section II G, we present some numerical estimates 



for the average || . ||t r distance of two randomly chosen density matrices. We first have to choose a measure over 
-B P os,i- All density matrices can be written as p = ^ \pu with Et=i = 1- The eigenvalues lie on a (N — 1)- 
dimensional simplex S in BL . We use the Euclidean metric || . H2 induced on the simplex. The Haar measure on the 
group of unitary matrices U(N) induces a uniform measure on the set of projectors {pii}^ =1 - Together this defines 
a measure Mb P03 1 @|- Within this measure, one can express the average distance between two density matrices pi 
and P2, using the unitary invariance of |j . \\t r , as 



[|| PI-P2 \\tr]M BposA = V ol(5)4(i7(JV)) J dU Jo dX ^ ■ ■ ■ d\ k ^ d^ . . . dp k 6(J2 t \ - l)£(E»Mi - !) 

Tr l Ej x jPjj ~ U Y , Hjpjjrt\. 

The values obtained by a numerical calculation of Eq. ( |A6| ) are tabulated in Table ||. 



(A6) 



APPENDIX B: PREPARATION OF THE BATH 



To prepare the state 



Pb,P = Pb,ft ® ■ • ■ ® (4,l3i ( B1 ) 



given H\, — Ei=i lif/2 <8 hi, we first calculate the eigenvalues and eigenvectors of each qubit Hamiltonian hi. We 
prepare the state 

nti(^ e -°|o)(o| + e ^ e -|i>(i|)/^. (B2) 

with {eo,i,ei.i} the eigenvalues of qubit Hamiltonian hi. This can be done by changing an initial state |0)(0| with 
probability e _/3ei l /Zj into state |1)(1| for each i. We then rotate each qubit to its eigenbasis {|6jo)i l^ii)} : 

®ti^=®*=i(|6io)(0| + |6* 1 )<l|)- (B3) 
In total we perform 2k elementary qubit operations plus some constant classical overhead. 
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H a 




S 


B 


dimension 


N = 2, ..,2 4 


if = 2^,..,2 b 


N 


K 


locality 


c 3 = 4 


c b = 2 


4 


4 


sampling scale a 


1 


f(n,k,c s ,c b ) 


1 


1 



TABLE I. Some settings in the numerical simulation. 
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FIG. 3. An example of an unsuccesful equilibration for n — 1, k — 3 and (3 = 3. 
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FIG. 5. Means and median for n=l (500 samples) 
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FIG. 6. Means for n=2 (200-500 samples) 
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FIG. 7. Means for n=3 (50-100 samples). 



33 



0.4 



0.3 



S 0.2 - 



£, 0.1 



0.0 *- 



k=5 
k=6 



8 10 12 



0.08 

0.06 

:- 

Q 
CD 

ra 0.04 -N. 

> 
< 

0.02 
0.00 



k=5 
k=6 



0.15 



0.10 



0.05 



0.00 



2 4 6 8 10 12 

P 

FIG. 8. Means for n=4 (15-20 samples). 




10 12 



34 



dim N 


mean 


standard errors ^/var/(n — 1), n = 1000 


4 


0.90388 


0.00740588 


8 


0.96190 


0.00514057 


16 


1.00294 


0.00341226 


32 


1.01452 


0.00220363 


64 


1.02617 


0.00132233 



TABLE II. The average distance between two randomly selected density matrices. 
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